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We present an extension of the time-dependent Density Matrix Renormalization Group (t- 
DMRG), also known as Time Evolving Block Decimation algorithm (TEBD), allowing for the com- 
putation of dynamically important excited states of one-dimensional many-body systems. We show 
its practical use for analyzing the dynamical properties and excitations of the Bose-Hubbard model 
describing ultracold atoms loaded in an optical lattice from a Bose-Einstein condensate. This allows 
for a deeper understanding of nonadiabaticity in experimental realizations of insulating phases. 
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I. INTRODUCTION 

Properties of many-body quantum systems are difficult 
to compute, because the eigenstates are usually unknown 
and the large size of the Hilbert space makes even ap- 
proximate descriptions difficult to manipulate. Some nu- 
merical methods - such as variational methods or mean- 
held theories - are intrinsically approximate yielding sim- 
ple expressions describing most of the interesting phys- 
ical properties, especially in the thermodynamic limit. 
The so-called "quasi-exact" methods do not have such 
restrictions, and are supposed to converge to the exact 
quantum result when some parameter tends to infinity. 
Those include exact diagonalization methods and quan- 
tum Monte-Carlo methods P, (the limiting parame- 
ter being the statistical sampling). For one-dimensional 
(ID) systems, particularly effective is the Density Ma- 
trix Renormalization Group (DMRG) approach 0, 0|. 
The states, in modern implementations of the method, 
are represented by the so called Matrix Product States 
(MPS) (for a recent smooth introduction, see Q). In nu- 
merous quantum systems, the ground state is very well 
approximated by a MPS, allowing its efficient descrip- 
tion with a relatively small number of parameters, which, 
moreover, does not increase exponentially with the sys- 
tem size. To some extent, DMRG shares the advantages 
of both variational and quasi-exact methods. 

While the original application of these methods con- 
sidered stationary properties, recent years brought a big 
progress in addressing the dynamics. Arguably, it was 
the formulation of the Time Evolving Block Decimation 
(TEBD) algorithm enabling study of real time dy- 
namics of ID systems of a reasonable size, which triggered 
a rapid progress in the field. Soon, a time-dependent 
DMRG (t-DMRG) method was proposed @, H with a 
fast realization of a close equivalence between TEBD and 
t-DMRG. Various applications of the method addressed 
dynamics of quantum phase transitions (9l-[rT|. quantum 
quenches and thermalizations [l^ - [l5| . periodic driving 



[l6l. \vf\, to give just a few examples. 

As soon as the evolution of the system ceases to be 
adiabatic, the excited states become important in the dy- 
namics. The global properties of whole groups of excited 
states may be obtained from evaluation of e.g. the spec- 
tral function. Following the pioneering work of Hallberg 
[l8[, a lot of effort has been put in developing efficient 
schemes for spectral functions evaluations (see e.g. (l9l - 
[23| and references therein). 

In the present paper, our aims are somehow comple- 
mentary. We show how one may obtain selected, dy- 
namically important excited eigenstates of strongly cor- 
related systems in a relatively simple way. Instead of 
addressing spin chain systems, as typical for tests of var- 
ious DMRG improvements, we shall immediately turn 
our attention to cold atoms settings and consider a Bose- 
Hubbard model as realized by cold atoms in ID optical 
lattice potentials. In (ll| . we have shown on such a spe- 
cific example how to extract the energies of the few lowest 
excited states of the system, with, however, no access to 
their properties. In the present paper, we explain how 
to use the TEBD (or t-DMRG) method in order to ex- 
tract some of the excited states with quite good accuracy. 
Those states are expressed also in terms of MPS making 
the calculations of expectations values and correlation 
functions easy. While the area law theorem [24j assures 
that the ground state may be efficiently represented in 
MPS representation this is by no means guaranteed for 
a generic excited state. Our method will work only for 
those states for which such an efficient MPS representa- 
tion exists. We shall show, however, that this restriction 
leaves plenty of interesting ground to explore as, most 
often, many excited states up to a reasonably large ex- 
citation energy, are also efficiently expressed as MPSs. 
Interestingly, the method selects states that are impor- 
tant in the systems' dynamics. It seems more suited fol- 
low lying excited states (although we give other exam- 
ples too). Still, in many-body systems, the density of 
states increases usually very fast with energy, and even 
relatively small systems (such as e.g. encountered in ul- 
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tracold atomic gases) have plenty of excited states lying 
not far above the ground state energy. 

While we use as an example the Bose-Hubbard (BH) 
model describing an ultra-cold atomic bosonic gas loaded 
in a one-dimensional optical lattice, the method proposed 
is quite general and can be used for any system where 
TEBD (t-DMRG) works for sufficiently long times. The 
latter restriction is necessary since the method is based 
on Fourier Transform techniques with resolution being 
dependent on the integration time. We shall show later 
that the method has an unexpected self purifying prop- 
erty allowing sometimes to use much longer integration 
times than anticipated in standard applications of the 
algorithm [H, 0]. Possible future improvements of the 
method are discussed in the concluding section. 

In the following sections, we first describe the method 
in some details. Then, we apply it for extracting se- 
lected excited states populated either by the periodic 
driving or by the non adiabatic dynamics in the disor- 
dered Bose-Hubbard system. The illustrative examples 
of excited states show the power and the limitations of 
our approach. 



paves the way to study the properties of individual ex- 
cited many-body states. 

To evaluate \(pr(Ei)), the integral is approximated by 
a discrete series. To sum the series, a procedure to per- 
form the sum of two MPSs as a MPS is required, which 
we now briefly sketch. Any state of a ID system of M 
sites may be expressed in the MPS representation 0, [25[ : 
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where r^'*' are matrices and A^ vectors and \ij) span 
a local Hilbert space on site j. When the T's and the 
A's satisfy simple orthogonality relations detailed in [2(| , 
their physical interpretation is simple. For example, the 
spectrum of the reduced density matrix associated with 
bipartite splitting L : R = {1, . . . , 1} : {I + 1, . . . , M}, 

Tr R \tp)(ip\ is nothing but for a ; = 1,2.... The 

corresponding entanglement entropy given by 



(4) 



II. EXTRACTION OF EIGENSTATES BY 
FOURIER TRANSFORM 

The basic idea is to start from an initial state which is 
a "wavepacket", that is a superposition of various eigen- 
states, described by a MPS, to propagate it in real time 
with a time-independent Hamiltonian, finally to com- 
bine the MPS states at various times to reconstruct 
the excited states. The evolution of a state by a time- 
independent Hamiltonian with eigenbasis | e^) is given by 

|V>(*)) = 2Ci ex P(-*^A) c il e *)- where c t = (ei\4>(0)), 
The Fourier transform (FT) of the autocorrelation func- 
tion C(t) = M0)|V(t)) = EJc^expH^/ft) over 
a time interval T is the autocorrelation spectrum: 
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where sine a; = smx/x. In the limit of long time T, it 
yields narrow peaks at the Ei's with weights \ci\ 2 , that 
allows to determine the energy spectrum as shown in [ll| . 

Here we show that it is actually possible to extract 
from the dynamics also the excited eigenstates with large 
overlap - those that contribute most to the dynamical 
wavepacket. We perform a FT directly on the many-body 
wavefunction \4>(t)) : 
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For long T, \4>r(Ei)) — > Cj|ej), providing us with the 
targeted eigenstate. The method selects excited states 
relevant for the dynamics, although myriads of other 
states may be present in the same spectral region. This 



is an important tool to characterize the strongly corre- 
lated systems, e.g. topological phases (27| . The collection 
of A values forms the entanglement spectrum. 

To describe exactly a generic state in terms of a MPS, 
a large number (exponentially increasing with M) of a/ 
values is needed. However, quite often, physical states 
created in many experiments are only slightly entangled 
so that A^j =12 are rapidly decaying numbers, which 
allows for introduction of a rather small cutoff x m au 
sums above, resulting in tractable numerical computation 
Q. Our approach is limited to such states only. 

Following [H| one has that: M(T,X) + M(T',X') = 
M(T © T', A © A')- The result does not satisfy orthogo- 
nality relations that are vital for achieving efficient MPS 
algorithms. The algorithm to restore them (2|| ensures 
that the sum is stored in a memory-efficient way. Apply- 
ing the above two-step algorithm to perform elementary 
addition suffices to obtain \%px{Ei)). The technical details 
are presented in the Appendix. 

Our goal is to obtain an eigenstate as pure as possible, 
not merely to detect its presence. This requires long time 
integration which may affect the accuracy of the TEBD 
algorithm 0, 0|. The integration time may be signifi- 
cantly shortened using few simple ideas which we now 
describe. 

For simplicity let us restrict to a subspace spanned by 
two eigenstates = a\eo) + /3|ei) as a special case. 
Performing the Fourier Transform as in eq. gives: 

\MEo))=a\e )+p(T)\ ei ), (5) 

with /3(T) = sine AE = E x - Eq. Clearly 
/3(+oo) = 0, therefore large enough T guarantees conver- 
gence of \iPt(Eq)) to | eo) - Another possibility is to choose 
T such that f3(T) = 0, i.e. a multiple of . This ensures 
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that l^r^o)) = l e o) despite relatively short evolution 
time. When other eigenstates contribute significantly to 
the wavepacket, there is no possibility to choose such a 
time T that contributions from other eigenstates van- 
ish simultaneously. However, if for some state \ip), and 



its eigenstate |ej ) there exists an eigenstate 



with 



nonzero overlap on \ip), such that \Ei —Ei 1 \ <C \Ei —Ej\ 
if j 7^ io,ii, then choosing T as a multiple of is an 
optimal choice. 

If f3 3> a in eq. ([5]), a realization of (3(T) <C 1 is more 
difficult. One solution is to restart the Fourier Transform 
evaluation from a partially converged result. Let us first 
define a map fE,T{ip) — iPt(E). One may verify that for a 
two-level system, the n fold composition of /e,t satisfies: 

Ie.tW = fifETr 1 ^)) = «|e ) + /9 (sine— J |e x ) 

This requires performing evolution for time nT and re- 
sults in a power- n decay of (ei|^)/(eo|V ; )- This observa- 
tion remains true in the full Hilbert space with almost 
no modification. 

Another way to obtain the power-rt decay of other 
eigenstates is to use a proper window function (29| . In- 
stead of performing FT as in © , one calculates 
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with writ) = 0.5 (l — cos (^fr)) , called a Hahn window, 
being one possibility. This window gives a 1 /T 2 conver- 
gence (instead of 1/T without a window) of (3(T). 

Further possible improvements will be apparent when 
discussing specific examples below. 



III. BOSE-HUBBARD MODEL 

We now illustrate the method using the Bose-Hubbard 
(BH) model (30l . [3l] | as realized in a gas of ultracold atoms 
in an optical lattice. We are aware of limitations of the 
tight-binding Bose Hubbard description of the system 
for deep optical lattices (32l . [33[. The BH model pro- 
vides, however, a realistic and commonly used simplifi- 
cation of the problem. A pioneer experiment [3~H ob- 
served a quantum phase transition (QPT) from a super- 
fluid (SF) phase to a Mott insulator (MI), and stimulated 
research on interacting many-body systems with external 
perturbations, for a review see |35| . A similar QPT was 
observed in a ID realization of the BH model (36|. A 
87 Rb condensate was loaded in a deep two-dimensional 
"transverse" optical lattice potential realizing an array of 
one-dimensional atomic tubes. A shallower optical po- 
tential along the tubes coming from a pair of counter- 
propagating beams V/Er = ssin 2 (27rx/A) (in units of 
the recoil energy Eji = h 2 /2mjn > \ 2 ) was superimposed 
on the system. After preparing the system at different 



potential depth, the lattice was modulated periodically 
at different frequencies measuring the absorbed energy 
(for details see [36]). In another realization of this essen- 
tially ID experiment, the effects created by disorder were 
analyzed [33|. We shall consider both these situations in 
the following. 



A. Periodic lattice oscillations 

Let us consider first the excitations created by a pe- 
riodic modulation of the lattice. The problem has been 
addressed in (l6| where the energy gain due to such a 
modulation was studied using a ID Bose-Hubbard model 
in a parabolic trap and TEBD for the time evolution, 
obtaining results which were in a qualitative agreement 
with the experimental s pec tra (36| . Recently, a novel ex- 
tended study appeared [l7| which analyzes in detail dif- 
ferent possible excitations in the same model using both 
t-DMRG and a perturbative approach. This later calcu- 
lation considers smaller and simpler system so, for peda- 
gogical reasons, we start our analysis for cases discussed 

m na. 

The ID Bose-Hubbard Hamiltonian reads: 
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where M is the number of sites, bj (frt) is the destruction 

(creation) operator of an atom at site j, nj = bjbj the 
number operator. The tunneling amplitude J and the 
interaction energy U depend on the lattice depth s (see 
below), ej = c(j — ro) 2 is the energy offset of a given 
site, due to the harmonic trap with ro the position of its 
minimum. When the lattice depth is modulated, s(t) = 
So[l + dcos(ut)], the energy absorption is enhanced at 
frequencies leading to resonant excitation of some excited 
states. 

The time-dependence of s(t) maps onto the corre- 
sponding time dependence of the tunneling amplitude 
J(t) and interaction strength U(t) using, as in fltil. [l7j. 
the approximate formulae due to Zwerger (38|: 



J(t) 



-.s(t) 4 exp(- 



and 



U(t) = 4V2^s 1/4 (t) 
A 



,V2 
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where recoil energy units are used and s± is the depth of 
"transverse" potential creating the ID tubes. Following 
[l7| . we denote by Uq the value of U in the absence of 
modulation (<5 = 0) and by Jo the corresponding tunnel- 
ing amplitude. Also, we take s± = 30, a s = 5.45nm, 
A = 825nm and the curvature of the harmonic trap 
c = 0.0123 as in 
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Figure 1: (color online) Absorption spectrum obtained for a 
8% modulation (i.e. 5 = 0.08), after integration time t = 
100 (smooth red line) and t = 300 (filled circles connected 
by the black line). The data for t = 100 are multiplied by 
three. The higher resolution obtained for longer times allows 
partially to separate various components due to individual 
excited states. Circles correspond to numerical data, black 
line connects them to guide the eye. The inset shows the 
average occupation numbers of the initial ground state of the 
system. 



Consider first a system of N = 36 particles (cor- 
responding to Fig. 2 of pj|) in a deep lattice with 
so = 15. The corresponding characteristic interaction 
energy is Uo w 0.714, while the tunneling amplitude is 
J /Uo « 0.01. 

For any state in MPS representation - be it the ground 
state, an excited state or a wavepacket - is it easy to 
compute expectation values of simple operators, such as 
the average occupation number (nj) at site I, as well as 
its standard deviation An; = yj (nf) — (ni) 2 . The ground 
state of the system is obtained by the TEBD algorithm 
with imaginary time evolution. The average occupation 
number of sites, displayed in the inset of Fig. [TJ shows 
two Mott plateaus with single and double occupancy. By 
simulating a periodic driving of the lattice using a real 
time TEBD evolution, we find the absorption spectrum. 
We take the modulation to be 8%, that is 8 times bigger 
than in [l7|; this creates bigger excitation in the system, 
making extraction of the excited states easier. For the 
same purpose, we shift slightly the lattice and the bot- 
tom of the harmonic trap taking ro = 18.62345. This 
breaks the reflection symmetry of the system with re- 
spect to the center of the trap, making the analysis and 
visualization of the excitations easier. Still, the spectrum 
for shorter evolution time, t = 100 (for ft = 1, the unit 
of time is the inverse of the recoil energv) resembles the 
perturbative theory curve in Fig. 2 of [l7j. Longer inte- 
gration time t = 300 gives a better frequency resolution 
and several sharp peaks emerge. The absorption for time 
t = 100 multiplied by a factor of 3 seems to reproduce 
quite nicely the global amplitude of the broad structure 



around lo = Uo, showing that, for such times, absorption 
in this frequency region is linear. This is not true for 
isolated peaks at low frequency as well as for the peak 
around u = 1.75t/o, suggesting that these peaks corre- 
spond to isolated single levels (with the corresponding 
Rabi-type evolution). The integration is performed with 
time step At = 0.02 using a standard second order Trot- 
ter approach assuming maximal occupation at the sites 
imax — 6 and x = 50 (those values are highly sufficient 
for such a deep lattice and mean occupation of sites less 
than four). 

The absorption spectrum now serves as a guide for 
extracting the excited states. Suppose we are inter- 
ested in excited states contributing to the first small 
peak around oj/Uq = 0.24. We take the wavepacket 
obtained at t = 300 with the modulation at frequency 
uj/Uq = 0.24 as a starting wavepacket for the analysis. 
We first evolve this wavepacket for T = 10000 with time 
step At = 0.05 using a time independent Hamiltonian 
(i.e. <5 = 0) and find the autocorrelation function Ct(E). 
The latter shows a doublet around energy E/Uq ps 0.24, 
corresponding to two local excitations, one on the left side 
and one on the right side of the trap center. From this 
autocorrelation function, we obtain accurate eigenener- 
gies. In order to extract the eigenstate with the chosen 
eigenenergy E, we perform a Fourier transform on MPS 
states during a time sufficiently long for all states out- 
side the peak at E/Uo « 0.24 to be significantly damped 
(in this specific case, T\ = 3000 is convenient). This 
provides us with a wavepacket with a large overlap 
with the desired eigenstate (typically 0.7-0.9, while the 
overlap of the eigenstate with the initial wavepacket was 
in the range 0.0001-0.01). The last step is to "purify" 
the wavepacket \ip\) recursively: the wavepacket is 
propagated for a relatively long time (Tj = 6000 in our 
specific case). The Fourier transform of the time series 
at energy E again which yields a new wavepacket \ipi + i), 
where the weight of all other eigenstates has been de- 
creased (a Hahn window is used to speed up this de- 
crease). This leads to an exponentially fast convergence 
of the wavepacket \ipi) towards the desired eigenstate. We 
observed that numerical inaccuracies may cap the expo- 
nential convergence of this procedure after 15-30 steps. 
In all propagating steps, we have again used time step 
At = 0.05. In general, exact computational setup is ad- 
justed taking into account the initial overlap, details of 
the spectrum and numerical stability. 

Altogether, one iteration takes about 4 hours, using 
8 threads on 4 core Intel i7 processors (the program is 
multi-threaded). We have performed 6 iterations to check 
the accuracy; a single additional iteration suffices for a 
1% accuracy on the occupation numbers. 

As mentioned above, two closely spaced in energy 
states correspond to the absorption peak at oj/Uq = 0.24, 
they differ from the ground state by the hopping of one 
particle between two Mott domains with n — 1 and n = 2 
occupations occurring either on the left or on the right 
side of the trap center. Had the trap be placed symmet- 
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Figure 2: (color online) Average occupation numbers per site 
for some selected excited states of the system analyzed in 
figure Q] The black line with circles is the ground state, while 
red lines with squares show selected excited states. One of the 
two states corresponding to the absorption peak at uj/Uq = 
0.24 is shown in panel (a) (the other state is symmetric w.r.t. 
the trap center). Panel (b) gives the density for one of the 
states associated to the very small peak around ui/Uq = 0.37. 
Panels (c) and (d) present two states associated to the "high" 
frequency peak at uj/Uq = 1.75. 



rically with respect to the optical lattice, the two elemen- 
tary excitations would have been degenerate and the evo- 
lution from a symmetric ground state under a symmetric 
Hamiltonian would yield the symmetric combination of 
two excitations. A system with broken mirror symmetry 
allows for separation of these excitations, the left one is 
shown in Fig.[2Ji. Fig.[2j3 shows a similar excitation where 
a particle hops from an occupied to an empty side (this 
time shown for a "right" excitation). The corresponding 
energy is larger since such a situation occurs further from 
the center, where the harmonic trap is steeper. Two other 
excitations shown in Fig. [2}; and Fig. [2ji are the "left" and 
"right" excitations associated to the high frequency peak 
at uj/Uq = 1.75 and correspond to a particle hopping 
from single to doubly occupied Mott zone creating a site 
with triple occupancy. 

All these excitations were correctly identified in [l7|, 
based on energy considerations (possible for simple situ- 
ations occurring in the deep insulating regime for large 
lattice depths). Our excited states provide a direct veri- 
fication. Similarly, we have checked that the broad struc- 
ture around uj/Uq = 1 (in fact consisting of several dis- 
tinct peaks as visualized by t = 300 absorption spectra 
in Fig. [1]) corresponds to single particle-hole excitations 
in the central Mott domain as identified by [l7j . 

Fig. [3] presents another set of excited states. Fig. [3J: 
resembles a typical particle hole excitation in the central 
Mott domain, but not between neighboring sites. Thus, 
the corresponding excitation energy lies not in the broad 
structure but on its wings. All other plots in this figure 
present excited states not analyzed in (l7j . They corre- 
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Figure 3: (color online) Average occupation numbers for an- 
other set of excited states of the system analyzed in Fig. Q] 
The black line with circles represent the ground state while 
red lines with squares show selected excited states. 



spond to absorption peaks occurring on the left shoulder 
(a-c) and on the right (d-f) shoulder of the central struc- 
ture. They show either particle-hole excitations occur- 
ring in the lower n = 1 Mott domain or more complicated 
situations. 

It is apparent that our method provides reliably ex- 
cited states in the localized regime. All these excited 
states are well represented by an MPS with relatively 
small x- To discuss a more challenging situation, we 
move to another example, the absorption spectra stud- 
ied in jT(| where, due to a shallower optical lattice, some 
excitations may occur in the superfluid regions. In par- 
ticular, two structures in the absorption spectra around 
uj/Uq w 1 as well as uj/Uq ~ 2 — 2.3 were observed, in 
agreement with the experiment j36| . The former reso- 
nance corresponds to the gap in the MI phase (i.e. to the 
particle- hole excitation in the Mott phase). Its appear- 
ance was interpreted as an evidence for the presence of 
the Mott phase. The latter has no such easy explanation; 
it was argued heuristically [IB] that it corresponds to a 
particle- hole excitation in the superfluid component. 

Our method allows to test the correctness of this claim. 
We take parameters of Ref. [16] with a slight modifica- 
tion. As in the previous example, we break the symmetry 
of the problem by shifting the lattice minima with respect 
to the harmonic trap. Such a situation seems generic 
from experimental point of view and allows to separate 
individual excitations. The absorption spectrum, similar 
to that of Fig. 4 of Ref. [H], is shown in Fig. gj The 
initial state - ground state at Sq — 12 - is excited by a 
8 = 0.2 (20%) modulation for a duration t = 50/Er - 
the result is shown as a thick dashed line. We observe 
a broad peak around modulation frequency uj/Uq = 1 
as well as a much smaller broad peak around lu/Uq = 2 
[47j . Fig. 0] shows the absorption spectra for a more gen- 
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Figure 4: (color online) Absorption spectrum for a system of 
65 particles on 40 sites, in an harmonic trap with curvature 
c = 0.006, centered around ro = 20.62345. The mean lat- 
tice depth is so = 12 (corresponding to Uo = 0.6752 recoil 
energy). It is modulated sinusoidally with 1.67% relative am- 
plitude. The smooth black line with open circles correspond 
to a modulation duration t = 50 /Er; for longer modulation 
time t = 200 / Er, sharper peaks become resolved (red thicker 
line). The absorbed energy is divided by 4 in the later case, 
to compensate for a longer perturbation. The thick (green 
online) dashed line shows the absorption curve for a much 
stronger 20% modulation for t — 50/ Er, rescaled by the in- 
tensity factor 12 2 . 



tie driving, i.e. 1.67% modulation. For an easier com- 
parison, the strong modulation curve is rescaled by 12 2 
factor (square of the relative amplitude of modulation) . 
The comparison of the two curves show that the struc- 
ture around uj/Uq = 2 behaves perturbatively, while the 
broad resonance around uj/Uq = 1 shows signature of 
saturation. 

The wavepacket obtained for a given frequency, at the 
end of the modulation, is a starting wavepacket for our 
procedure. The strong 20% modulation, as used in [l6j . 
leads to a "too excited" wavepacket for the TEBD al- 
gorithm to work effectively over the long time required 
for a high-resolution FT (see, however, a discussion be- 
low) . Of course, the excited states do not depend on the 
strength of the modulation (while the shape of the ab- 
sorption curve in the nonperturbative regime does). For 
that reason, we use the wavepacket obtained after a "gen- 
tle" modulation with a smaller (1.67%) amplitude. 

Fig. [4] shows also the absorption curve for a longer 
modulation interval t = 200/ Er which allows for a par- 
tial resolution of broad peaks. In particular, the struc- 
ture around uj/Uq = 2 shows narrower partially resolved 
peaks around uj/Uq = 1.84, 2.06, 2.25. 

Let us concentrate first on the wavepacket modulated 
with frequency uj/Uq = 2.3. We perform a FT of the 
autocorrelation function over a long time and observe 
several partially overlapping peaks. As shown in the in- 
set of fig. [SI we identify a doublet of states, relatively 
well separated from the others, significantly excited from 




Figure 5: (color online) Black circles: average occupation 
number for the ground state of 65 particles in a lattice with 
depth a — 12 (40 sites) and an harmonic trap with curva- 
ture c = 0.006, centered around ro = 20.62345 (to separate 
left and right excited states w.r.t. to the trap center). The 
filled squares (red online) connected by a dashed line and the 
crosses (green online) connected by a dotted line show occu- 
pation numbers for two extracted excited states showing a 
particle-hole excitation in the SF region. The excess energies 
(energies measured with respect to the ground state energy) 
of the states are AE/U = 2.2502 (resp. AE/U = 2.2434) 
for the "red squares" (resp. "green crosses) in excellent agree- 
ment with the positions of two closely spaced peaks in the 
FT of the autocorrelation function shown in the inset. The 
lower panel shows the standard deviation of the occupation 
number on each site, for the ground state (open circles) and 
the excited states, the latter shows reduced fluctuations. 



the ground state by the modulation. Fig. [5J shows the 
average occupation numbers (ni) at site I for the two 
excited states, extracted by our method, together with 
the ground state. Observe that, indeed, these states cor- 
respond to a particle-hole excitation in the SF regime 
(with average occupation numbers between 2 and 3) con- 
firming the suggestion of [ll|- The two states are quasi- 
symmetric images; they have slightly different energies, 
because of the choice ro = 20.62345 which breaks parity 
w.r.t. the trap center. 

We have found that, in order to obtain a fully con- 
verged excited state (e.g. the excitation appearing on 
the right hand side of the center of the trap repre- 
sented as crosses in Fig. [5j, another technique is use- 
ful. Knowing from the preliminary runs where the lo- 
cal (in space) excitation should be created, we pre- 
pare a new wavepacket by a local modulation, i.e. lo- 
cal chemical potential is modulated for a few nearby 
sites (around site 25 for this particular case). Explic- 
itly Hi{t) = fii(l + <5exp[— (i — iq)/w 2 ] sin(wi)) where w 
is the spatial width of the modulation. Experimentally, 
such a modulation could be in principle created by an 
additional focused laser. Here it is used in the numerical 
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simulation to speed up the convergence as well as to iso- 
late excitations quasi-degenerate in energy but spatially 
well separated. 

The data in Fig. [S] show some interesting differences 
when compared to the Mott regime previously discussed. 
While, in Figs. [2fc and [2ji. the excited states display a 
transfer of particles towards the trap center - explaining 
why they have an excess energy (compared to the ground 
state) smaller than 2C/ - it is the opposite for states in 
Fig. 03 where particles are transferred further from the 
trap center, and the associated excess energy is larger 
than 2Uo. The MPS representation of the excited states 
makes it easy to calculate several observables, not only 
the average occupation number of sites. The lower plot in 
FigJS] shows the standard deviation of the occupation of 
sites An; , both for the ground state and the excited one. 
A small (resp. large) standard deviation is a signature of 
a Mott insulator (resp. superfluid region). The excited 
state has lower An; in the region of excitation, proving 
that it is more an "excited Mott". This is confirmed by 
the occupation numbers themselves which are close to 
integer values for the excited state, with larger fractional 
part for the ground state. 

To explore other excited states contributing to the 
structure around lo/Uq = 2, we prepare the wavepacket 
with modulation frequency lo/Uq = 2. The first attempt 
to calculate the autocorrelation function by propagating 
such an excited wavepacket with x = 70 (and a local 
Hilbert space on each site containing occupation num- 
bers up to imax=7) failed, as e.g. the energy in the 
system was not conserved. We chose, however, to simply 
run the algorithm for a much longer time and observed 
an unexpected behavior depicted in Fig. [5] 

Fig. [5k presents the energy of the wavepacket versus 
time (in recoil units, multiplication by Uq = 0.6752 give 
time in units of 1/Uq). After the initial decrease, we ob- 
serve stabilization of the energy at some value above the 
ground state. Fig. presents the corresponding time 
dependence of the entropy of entanglement in the system 
defined as S = max; 5; with 5; given by Eq. It also 
decays and stabilizes. This behavior can be understood 
by performing the FT of the wavepacket autocorrelation 
function (shown in Fig. |Bt). While the FT over the time 
interval [0,5000] displays broad erratic structures, the FT 
over the time interval [15000,20000] gives sharp peaks at 
some frequencies. The heuristic understanding of this be- 
havior is the following. The initial excited wavepacket is a 
linear combination of eigenstates of the system. Some of 
them are well represented as an MPS, some other ones are 
represented poorly. During the time evolution the "well 
represented" part evolves accurately while the remaining 
part produces signatures of TEBD algorithm breakdown. 
During the course of evolution, a self-purification of the 
wavepacket occurs with the non-converged fraction seem- 
ingly disappearing. 

The information on the non-converged fraction is ac- 
tually lost in the TEBD evolution, when some Aa val- 
ues are discarded (those with a > X-) No, or little, in- 
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Figure 6: (color online) Long time evolution of a wavepacket 
prepared by a modulation at frequency uj/Uo = 2 and then 
evolved with the TEBD algorithm. The excess energy over 
the ground state is plotted in top frame (a). It is not con- 
served but, after the initial drop, it stabilizes. Similar behav- 
ior occurs for the entanglement entropy (b), with additional 
oscillations because the wavepacket is a sum of relatively few 
localized states. The FT of the autocorrelation function (c) 
does not display any clear peaks if performed in the [0,5000] 
time interval (black thin line), but definitely shows several 
well defined excited states over the time interval [15000,20000] 
(thick line, red online). 



formation is lost on the well represented excited states. 
Thus, unconverged TEBD time evolution filters out in- 
formation which anyway cannot be obtained, purifying 
the wavepacket. 

A further confirmation comes from the fact that the 
peaks of the FT provide the energies of excited states 
which can now be extracted. FigJTJshow examples of ex- 
cited states involving excitations in the superfluid regime. 
Contrary to the previous examples, these excited states 
show increased particle number fluctuations, see lower 
panel of FigfT] 

States presented in Fig. [5^-Fig. [5J: are other exam- 
ples, well localized on the edge between the (n;) = 1 and 
(n;) = 2 zones. The corresponding standard deviations 
shown in Fig.[H]show that these excitations, nevertheless, 
lead to an increase of local particle number fluctuations. 
Fig. [8]i shows a state with much larger excitation en- 
ergy, almost &Uq. It shows that our method may produce 
states with quite high and complicated excitations. This 
state displays several particle transfers as well as almost 
a transfer of a part of the right hand side Mott (n;) =2 
plateau into a superfuid region. 

In effect a vast majority of states identified by the FT 
of the autocorrelation function can be extracted. Natu- 
rally, most of these states correspond to local excitations 
and not too high entanglement entropy. 
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Figure 7: (color online) Average occupation numbers of exem- 
plary excited states showing a significant superfluid-type ex- 
citation, when compared to the ground state (black open cir- 
cles). Red squares connected by the dashed line correspond to 
the excited state with excess energy AE/Uo = 2.2418 (lower 
energy peak of the doublet in the inset of Fig. [EJ. Green 
crosses connected by dotted lines correspond to a state (excess 
energy AE/Uo = 1.9504) with significant excitation around 
the trap center. The lower panel shows the standard deviation 
of the occupation of sites for the ground state (open circles) 
and two excited states. The latter show regions of increased 
superfiuid-like fluctuations. 




Figure 8: (color online) Average occupation numbers of ex- 
emplary excited states showing localized particle-hole exci- 
tations, compared to the ground state (black open circles). 
Red squares connected by dashed lines correspond to the 
excited states. (a) Excitation similar to those discussed 
for the smaller system, i.e. particle-hole excitation on the 
transition region between two Mott plateaus (excess energy 
AE/Uo = 1.8417). (b) and (c) show two particles being trans- 
fered between sites, the excess energies are AE/Uo = 1.9988 
and AE/Uo = 2.0539. (d) gives an example of higher lying 
excited state AE/Uo = 7.5280 with a partially melted, right 
(n{) —2 Mott plateau. 




Figure 9: (color online) Standard deviations of the occupa- 
tions number, for the states presented in Fig. [8] Local exci- 
tations lead to an increase of particle number fluctuations. 



B. Florence experiment 

A second important example is the ID experiment of 
the Florence group [37} where a second, weaker lattice 
potential (incommensurate with the primary lattice) is 
added to simulate a disorder. Its strength is denoted 
by S2- Theoretically, the addition of a disordered poten- 
tial reveals a new insulating, compressible and gapless 
phase, called Bose glass (BG) (3lj . Experimental results 
confirm its existence [13, HH, but are not fully conclusive: 
the system is prepared from a low temperature BEC by 
ramping up an optical lattice, starting from a SF ini- 
tial state. Preparing a MI or a BG requires, generally 
non-adiabatic, passage through a quantum phase transi- 
tion [ZoL I4l1 | , thereby exciting the system. 

The Florence experiment is described in [3^ . We simu- 
lated the experiment as detailed in (Tl| : we integrate nu- 
merically the evolution of the system when s increases us- 
ing the TEBD algorithm Q , starting from the SF ground 
state at s = 4. At the same time, a secondary shallow 
lattice creating the disorder is turned on. Its effect is to 
modify the energy offset of the harmonic wells at different 
sites yielding now: 



ej = c(j - r ) 2 + s 2 E R2 sin 2 



TTjA 

A 2 



(10) 



The curvature of the trap, c depends on the harmonic 
trap and on the additional trapping potential created by 
the transverse profile of the laser beams [37| . 

As s increases, we found [ll| that the overlap of the 
dynamically evolved wave packet on the ground state for 
that s, rapidly decreases in the region of the SF-MI tran- 
sition (around s = 8 — 9 for the experimental parameters) 
down to less than 10% for the final s = 14, a signature of 
broken adiabaticity. The situation is even worse (with the 
final overlap of the order of 10~ 6 ) for the disordered sys- 
tem. Disorder enriches the energy level structure around 
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Figure 10: Autocorrelation spectra, Eq. fT}, obtained dy- 
namically for s=14 after switching on the optical lattice, for 
disorder strengths s 2 = (a) 0, (b) 0.4375, (c) 2.1875. All pa- 
rameters are taken to approximate the experimental situation 
[37l | . The energy levels of the system appear as peaks (origin 
at the zero-point energy), with height \ci\ 2 from eigenbasis 
expansion. Peak labels are for further reference. 



the ground state, killing adiabaticity both for shallow [42| 
and deep [ll| lattices. 

To identify the induced excitations, we compute the 
various excited states building the dynamically created 
wavepacket, by evolving it further at the constant final s 
value and calculating the FT in Eqs. (|1I2[) , as explained 
in the previous Section. In particular we characterize the 
states obtained using similar observables (easily obtained 
in the MPS representation [28|]): the average occupation 
number (n;) at site / and its standard deviation An;. As 
a by-product, we can also read off the average number of 
particles in the left part of the system TV; = 5_3<=i ; n i 
and its standard deviation ANi. This quantity is rela- 
tively easy to measure in a experiment. In our system 
it qualitatively resembles the entanglement entropy, as 
shown below [43|. 

Fig. [PHI shows the autocorrelation spectrum, eq. JT]), of 
the dynamically created wavepackets, at increasing disor- 
der strengths. In the absence of disorder (S2 = 0), about 
ten states are significantly excited. In Fig. [TTJ we show 
various relevant quantities for the ground state, few ex- 
cited states and the wavepacket. The average occupation 
number (m) on each site I has the well known "wedding 
cake" structure, with large MI regions with integer (n;) 
separated by narrow SF regions. Because the energy ex- 
cess brought by non adiabatic preparation is small, all 
significantly populated excited states have similar shapes. 
Clearly, all excitations take place in or around the SF 
regions: these are transfers of one atom from the edge 
of a Mott plateau to the edge of another Mott plateau 
or to the neighboring SF region ("melting" of the Mott 
plateau). 

The standard excitations in an homogeneous system 
(without trapping potential) such as a particle-hole ex- 
citation in a Mott plateau are absent in the stationary 




Figure 11: (color online) Properties of states of the trapped 
BEC in a deep optical lattice (s=14), without disorder (s2=0), 
ro=0. 12345 in Eq. (|10|) to break the parity symmetry w.r.t. 
the trap center. The black line refers to the ground state, the 
red line to the dynamically prepared wavepacket, the brown 
(with crosses) and blue lines to the two excited states with 
the largest populations denoted as 1 and 2 in Fig. 110b . (a): 
Average occupation number (ni) on each site, (b): Differ- 
ence Dl = (m) - (n;) G nd state (c) : All; = y/ {nf) - (n t ) 2 . 

(d): Entanglement entropy, Eq. (Q; almost zero in the Mott 
plateaus for eigenstates - implying their approximate separa- 
bility; large for the wavepacket. (e): Variance of the number 
of atoms to the left of site I. It is akin to the entanglement 
entropy, showing a marked difference between the eigenstates 
and the wavepacket. 



states because they are dynamically unstable. They are 
also absent in the wavepacket because they are energet- 
ically too costly: a particle-hole excitation costs at least 
one interaction energy U ~ 0.6-Er, much larger than the 
excess energy of the wavepacket 0.1-Er. Observe, there- 
fore, that we are in a regime significantly different from 
the one discussed in the previous Section where excited 
states populated by modulation of the lattice depth have 
been considered. Here we are much "closer" in energy to 
the ground state. 

The description used in [44| where the ground state 
is contaminated by local particle-hole excitations is not 
compatible with our findings. A description in terms of 
melting of the MI (4f| seems more relevant. The local oc- 
cupation fluctuation An; confirms the existence of large 
MI regions with low An; separated by SF peaks with 
larger An;. While the ground state and excited states 
look very similar, the wavepacket is different, with higher 
An; in the SF peaks. An even greater difference is re- 
vealed by the entanglement entropy. For the eigenstates, 
it is essentially zero in the MI regions and displays sharp 
peaks in the SF regions. In stark contrast, the Si of the 
wavepacket is also non-zero in the Mott plateaus. Sim- 
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Figure 12: (color online) Same as Fig. QT] but for small dis- 
order (s 2 =0.4375, r =0 and 0=0.5432 in Eq.(JlO|). Brown 
(blue) lines correspond to excitations A and B indicated in, 
Fig. HOf b). Mott insulating regions are clearly visible for all 
eigenstates, more vague for the wavepacket. 



ilarly, ANi is much larger in the Mott plateaus for the 
wavepacket than for the eigenstates. This suggests the 
existence of long-range entanglement between SF inter- 
faces linking different MI regions caused by long range 
particle dislocations, absent in the eigenstates. 

Disorder strongly modifies the properties of the sys- 
tem. Possible phases are analyzed in details in (4(|. We 
consider first a small disorder ,S2=0.4375. The breakdown 
of adiabaticity is stronger in presence of disorder, with 
twice larger excess energy - still below the particle-hole 
excitation energy - and more states significantly excited. 
The properties of various states are shown for small S2 
in Fig. [T2l Several MI regions - identified by plateaus in 
(ni) - still exist, separated by intermediate parts which 
are either SF or a BG. Remarkably, the ground state and 
the excited states do have very similar structures, with 
several visible Mott plateaus. In contrast, these plateaus 
are less visible (only the central one seems to survive, 
with a reduced size) for the wavepacket. This can be in- 
terpreted as a partial melting of the BG and MI phases, 
producing a thermal insulator (45| . 

For strong disorder, = 2.1875, no MI state exists 
and the ground state of the system forms a BG. The 
exemplary excitations are shown in Fig. 1131 The occupa- 
tions of various sites strongly fluctuate, low lying excita- 
tions seem quite similar to the ground state. Excitations 
are local in character modifying the occupation number 
and its variance in selected (not always close) sites only. 
The entanglement entropy is small, except in small SF 
pockets. Again, the wavepacket has large entanglement 
entropy and An; , with numerous peaks indicating melted 
regions. 




20 40 , 60 80 100 



Figure 13: (color online) Properties of states of bosons 
trapped in a deep optical lattice (s=14) in the presence of 
strong disorder s 2 =2.1875. Here r =0.1075 and 0=0.0304 in 
Eq.([T0}). Brown (blue) lines correspond to excitations X and 
Y in Fig. HOf cl. The entanglement entropy is much larger 
for the dynamically created wavepacket than for stationary 
states, and An; has many more peaks, indicating a signifi- 
cant melting of the Bose glass. 



IV. CONCLUSIONS 

In conclusion, we have built an extension of the TEBD 
algorithm capable of extracting, from a non stationary 
state (created dynamically), the eigenstates dynamically 
populated, and the nature of the "defects" created. In 
the specific case of the BH model realized in a ultracold 
atomic gas, a quasi-adiabatic quench populates excited 
states that are essentially similar to the ground state. In 
contrast, the dynamical wavepacket has a markedly dif- 
ferent character, with a significant entanglement across 
the whole sample. By modulation spectroscopy, higher 
excited states may also be selectively reached and ana- 
lyzed. 

It is clear that our method will not provide an ac- 
cess to all excited states of a given many body system. 
By construction, only those excited states which are well 
represented by a MPS within a given assumed size of 
the internal Hilbert space x m &y be obtained. On the 
other hand, the TEBD algorithm seems to create a "self- 
purification" mechanism which enables to reach integra- 
tion times much longer than originally expected 0, 0|- 
Naturally the result of this integration is not a good 
approximation of the true dynamics but contains infor- 
mation about the well represented excited states. Thus 
our method seems to have larger range of applicability 
than expected on the basis of standard limitations of the 
TEBD algorithm. 

In that context, wavepackets with low excess energy 
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w.r.t. the ground state created via quasi-adiabatic 
quench simulating the Florence experiment with disor- 
der [53] behave exceptionally well. Higher excitations 
reached in modulation spectroscopy were more computer 
intensive and revealed that the method cannot be used 
as a kind of "black-box" approach but has to be tuned 
to individual cases. While this is certainly some draw- 
back, the method potential capabilities (as exemplified 
with the highly excited state shown in fig. EH) seem quite 
encouraging. 
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VI. APPENDIX: IMPLEMENTATION OF 
MPS- ADDITION 

The addition of two MPS, an elementary step in the 
FT, produces A vectors with increased size. Thus it seems 
that the x value required for representing successive sums 
should increase fast, ruining the efficiency of the MPS 
representation. This is not the case: the A^j =1 2 for the 
partial sums decrease relatively fast so that the small- 
est components can be safely discarded [2^|, preventing 
the cutoff to blow up, a key point for the success of 
the method. Altogether, several thousands states can be 
summed, at the price of roughly doubling the x value to 
preserve the accuracy of the MPS representation: when 
the series is dynamically created, as in eq. ([B]), the result 
is close to an eigenstate, i.e. simpler than the summands 
which are wavepackets, and the vectors added together 
are closely related. 

To be able to compute the sum of several thousands 
of MPS-vectors two issues have to be resolved: restor- 
ing the canonical form of a MPS vector and performing 
truncation in \- The former has been described in (26| . 
The truncation from x = Xi to x — Xi can be- done in 
two nonequivalent ways: "blunt" and "canonical". The 
canonical truncation is done in parallel with restoring 
the canonical form. The algorithm described in (26| pro- 



cesses sites sequentially restoring orthogonality relations 
at each of them one at a time. The "left" and "right" 
Schmidt vectors for site k are: 

\^ ak k ) =£ rfif AWrL 2 ].^ . . . r^% k \h, . . . , i k ) (11) 

ai °>k 

*1 *k 

\,k+l...M\_ r [fc+l],i fc+ i A [fc+l] v [M].iM I- • \ 

a fc+1 ,...,a M 

(12) 

W=E A aifc fc )®l^ 1 '-' W > ( 13 ) 
<*k 

For a MPS vector in its canonical form, these vectors 
are orthonormal eigenvectors of the reduced density ma- 
trices pi 1 -'*], p[k+i-;M]_ xhe algorithm in [H] restores 
orthonormality of left and right Schmidt vectors at site k 
assuming it has already been applied for sites 1, . . . , k— 1. 
After the orthogonality has been restored at site k, and 
the A^ sorted by decreasing magnitude, we zero all \\% k 
for xi < a k < X2- The truncation error is estimated 

by Ea 2 fc=Xl +i( A S) 2 El- After all sites have been pro- 
cessed, at each site left Schmidt vectors are orthonormal 
eigenvectors of the appropriate reduced density matrix. 
Orthonormality of right Schmidt vectors at site k is lost 
when truncation of A^ +1 ] takes place. If the algorithm 
were applied again (without truncation), both left and 
right Schmidt vectors would be canonical for any site. 
We do not perform this second part since we do not use 
orthonormality of right Schmidt vectors for the addition 
procedure. 

The blunt way to truncate the MPS state is to first 
use the previous algorithm with no parallel truncation at 
all, only then perform truncation of the vectors \*- k ' at all 
sites. This method is mathematically flawed as, starting 
from the second site truncation, the canonical form is 
lost. However, when the x value used for FT evaluation 
is larger than x for the temporal evolution, both methods 
agree well. 

In that respect, wavepackets with low excess energy 
w.r.t. the ground state created via quasi-adiabatic 
quench simulating the Florence experiment with disorder 
|37j behaved exceptionally well. The converged results 
could be obtained with % = 80 as tested and compared 
with x = 120 or x = 150. 

In contrast, extraction of selected eigenstates from 
wavepackets created by harmonic modulation is more 
demanding. Too strong modulation creates wavepack- 
ets containing apparently strongly entangled contribu- 
tions. In such a case, a long time propagation of such a 
wavepacket is not possible even for large x values. Low- 
ering the amplitude of the modulation and adjusting its 
frequency is a key ingredient for successful extraction of 
a given excited state. The results presented in this pa- 
per have been obtained for x — 120 again testing the 
convergence by comparison with x = 150. 
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The convergence may be also tested by monitoring the 
growth of the smallest values. 
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